#!/bin/gawk -f
# qmodr: Script to determine contents of a restraint file.
# The script writes a modified restraint file, 
# inserting data identifiers (if not <ide>), 
# and comment delimiters (if not <cmm>).
# Input variables: see usage statement
# Exit code: 0 (ok), or 1 (error)
# If <ide>, <cmm> and <mr> are all 0 the script just copies the input file.
# =============================================================
# Ton Rullmann, October 8, 1996
# =============================================================
# Comment on the code:
# Data types are guessed from (1) the data (2) identifiers (3) comments
# (from high to low priority), unless <ide> equals 0. 
# Comments are read from all lines, including data lines (see guesscmm; 
# interpretation switches are set by <how>).
# Types derived from (2) and (3) are stored in <deftyp>, and their origin
# is stored in <deffrom>.
# Types derived from (1), i.e. the data itself, override <deftyp>. These
# types are temporarily stored in <typ> (see readatXYZ and guessdat).
# If <!nosavtyp> the <typ> value is saved in <deftyp>. 
# When SS or HB restraints appear amidst NOE restraints, without extra comment
# lines being present, the type can be switched back to NOE by setting
# <nosavtyp> to 1, or by setting <how> to 2 or 3. These options are not
# exactly equivalent; try which one works best.
# The type expressed on the last printed START_OF... identifier is stored
# in <dattyp> (see prtide).
# If an indication about upper or lower bounds is found in the comments
# (see guesscmm) a ...BOUND_DATA identifier is printed. The type is
# stored in <defuplo>.

BEGIN {
    logfil = "qmodr.log"
    if ( ARGC <= 1 ) {
	print "usage: qmodr filtyp=1|2|3|4|99 [mr=1] [ide=1] [cmm=1] [fatal=1] [nosavtyp=1] [how=value] file"
	print "   filtyp:   file type number (required)"
	print "              1 = DISGEO"
	print "              2 = XPLOR"
	print "              3 = BIOSYM"
	print "              4 = DISMAN/DIANA"
	print "             99 = user defined"
	print "   mr:       if 0, don't look for coordinate section of an mr file"
	print "             else, chop off all but first model and (if ide) insert identifier"
	print "   ide:      if 0, don't insert identifiers (i.e. assumed to be present)"
	print "             else, insert"
	print "   cmm:      if 0, don't insert comment delimiters (i.e. assumed to be present)"
	print "             else, insert"
	print "   fatal:    if 0, don't exit at conflicts of data type"
	print "   nosavtyp: if 1, don't save a type derived from data as new default"
	print "   how:      how to interpret comments or data:"
	print "             1: string DISTANCE -> interpret as NOE type"
	print "             2: set default type derived from comments on data line to NOE"
        print "               (this will annihilate effect from identifiers and comment lines)"
	print "             3: set default type derived from data to be NOE"
	print ""
	print "Error/warning output (if any) appears on " logfil
	exit ( 1 )
    }
    if ( !ide && !cmm && !mr ) {
	system( "cat " FILENAME )
	exit
    }
    cmmorig = cmm
    filtyp_max = 4
    filtyp_unk = 99
    if ( filtyp < 1 || filtyp > filtyp_max ) filtyp = filtyp_unk
    typ = "none"
    deftyp = "none"
    dattyp = "none"
    deffrom = "none"
    defuplo = ""
    dismanreslin = ""
    lastdata = 0
}

function countdatafields() {
    for ( i=1; i<=NF; i++ ) {
	if ( match( $i, /^[!#]/ ) )
	    return ( i-1 )
    }
    return ( NF )
}

function garbage () {
# store the line, and optionally insert comment delimiter
    if ( cmm ) {
	cmmsto[++cmmcnt] = "# " org
	print "* Made into comment: " org >> logfil
    }
    else
	cmmsto[++cmmcnt] = org
}

function guesscmm() {
    if ( !ide ) return
# Biosym: interpret only the Biosym identifiers
    if ( filtyp == 3 ) {
	if ( $0 ~ /^#NOE_DISTANCE/ ) {
	    deftyp = "NOE"
	    deffrom = "BIOSYM identifier"
	    return
	}
	if ( $0 ~ /^#DISTANCE/ ) {
	    deftyp = "DIS"
	    deffrom = "BIOSYM identifier"
	    return
	}
	if ( $0 ~ /^#NMR_DIHEDRAL/ ) {
	    deftyp = "TOR"
	    deffrom = "BIOSYM identifier"
	    return
	}
	if ( $0 ~ /^#/ ) {
	    deftyp = "OTH"
	    deffrom = "BIOSYM identifier"
	    return
	}
	return
    }
# Info whether we have upper- or lower bounds?
    if ( filtyp >= 4 ) { 
	df = ""
	if ( $0 ~ /UPPER/ ) {
	    df = "UPPER"
	}
	if ( $0 ~ /LOWER/ ) {
	    if ( df == "" )
		df = "LOWER"
	    else
		df = "UPPLOW"
	}
	if ( df != "" ) defuplo = df
    }
# How==2 option: reset deftyp
    if ( how == 2 ) {
	deftyp = ""
    }
# Now check for clues
    if ( $0 ~ /DIHEDRAL/ || $0 ~ /TORSION/ || $0 ~ /ANGLE/ || \
	 $0 ~ /[#! ]*PHI( |$)/ || $0 ~ /[#! ]*PSI( |$)/ || $0 ~ /[#! ]*OMEGA( |$)/  ) {
	deftyp = "TOR"
	deffrom = "comment"
	return
    }
    if ( $0 ~ /H[- ]*BOND/ || $0 ~ /HYDROGEN BOND/ ) {
	deftyp = "HB"
	deffrom = "comment"
	return
    }
    if ( $0 ~ /S[- ]*S[- ]*BRIDGE/ || $0 ~ /S[- ]*S[- ]*BOND/ || $0 ~ /S[- ]*S[- ]*(CON|RE)STR/ || $0 ~ /DI[- ]*SUL(F|PH)IDE/ ) {
	deftyp = "SS"
	deffrom = "comment"
	return
    }
    if ( $0 ~ /[#! ]*NOES?( |$)/ || $0 ~ /[#! ]*NOE-/ || $0 ~ /[#! ]*NOE_/ || $0 ~ /PROTON *- *PROTON/ || $0 ~ /INTERPROTON/ ) {
	deftyp = "NOE"
	deffrom = "comment"
	return
    }
#    if ( how == 1 && $0 ~ /DISTANCE *(CON|RE)STRAINT/ ) {
    if ( how == 1 && $0 ~ /DISTANCE/ ) {
	deftyp = "NOE"
	deffrom = "comment"
	return
    }
    if ( $0 ~ /VOLUME/ || $0 ~ /INTENSITIES/ ) {
	deftyp = "OTH"
	deffrom = "comment"
	return
    }
#    if ( $0 ~ /EXCHANGE *TIME/ || $0 ~ / HZ / ) {
#	deftyp = "OTH"
#	deffrom = "comment"
#	return
#    }
# How==2 option: default deftyp
    if ( how == 2 ) {
	deftyp = "NOE"
	deffrom = "comment"
	return
    }
}

function readatXPLOR() {
    if ( !ide ) return
# find atoms in specific fields - cf. convXPLORxxxrestr.old
    if ( filtyp == 2 ) {
	gsub( /\(/, " " )
	gsub( /\)/, " " )
	AA = $6
	BB = $11
# this does not catch all cases of course....
	if ( match ( $0, /SEGI/ ) ) BB = $14
# first line of torsion angle restraints split over two lines
	if ( NF == 11 ) typ = "TOR"
    }
}

function readatBIOSYM() {
    if ( !ide ) return
    sub( /^.*:/, "", $1 )
    AA = $1
    sub( /^.*:/, "", $2 )
    BB = $2
}

function readatDISMAN() {
    if ( !ide ) return
    if ( nfield < 7 ) {
	AA = $1
	BB = $4
    }
    else {
	AA = $3
	BB = $6
    }
}

function guessdat() {
# routine to guess from the data AND print the identifier
    if ( !ide ) return
# remove on-line comment (if present)  -- other types than XPLOR may use # ... problem
# not needed anymore ???
    if ( $0 ~ /!/ ) {
	sub ( /!.*$/, "" )
    }
# derive type from the atom names
# note: typ may have been set already in readatXYZ
# typ and atom names are initialized after reading a new line
# note: this might produce unwanted results when the data are for a 
# dihedral angle, but that is an unlikely event
    if ( ( AA ~ /(^N)|(^HN$)/ && BB ~ /^O/ ) || \
	 ( BB ~ /(^N)|(^HN$)/ && AA ~ /^O/ ) )
	typ = "HB"
    else if ( ( AA ~ /(^CB$)|(^SG$)/ && BB ~ /^SG$/ ) || \
	 ( BB ~ /(^CB$)|(^SG$)/ && AA ~ /^SG$/ ) )
	typ = "SS"
    else if ( ( AA ~ /^HB[123]$/ && BB ~ /^SG$/ ) || \
	 ( BB ~ /^HB[123]$/ && AA ~ /^SG$/ ) )
	typ = "SS"
    else {
	if ( how == 3 ) typ = "NOE"
    }
# now print the identifier (if necessary)
    prtide()
}

function guessdatUNK() {
# routine to guess from the data AND print the identifier - for UNKNOWN data type
    if ( !ide ) return
# remove on-line comment (if present)  -- other types than XPLOR may use # ... problem
# not needed anymore ???
    if ( $0 ~ /!/ ) {
	sub ( /!.*$/, "" )
    }
# derive type from the atom names (not fool proof - wants names isolated by spaces)
# note: typ may have been set already in readatXYZ
# typ and atom names are initialized after reading a new line
# note: this might produce unwanted results when the data are for a 
# dihedral angle, but that is an unlikely event
    if ( ( $0 ~ / H?NH?[0-9]* .* O / ) || \
	 ( $0 ~ / O .* H?NH?[0-9]* / ) )
	typ = "HB"
    else if ( ( $0 ~ / SG .* SG / ) || \
	      ( $0 ~ / SG .* CB / ) || \
              ( $0 ~ / CB .* SG / ) )
	typ = "SS"
    else if ( ( $0 ~ / HB[123] .* SG / ) || \
              ( $0 ~ / SG .* HB[123] / ) )
	typ = "SS"
    else {
	if ( how == 3 ) typ = "NOE"
    }
# now print the identifier (if necessary)
    prtide()
}

function prtide() {
# routine to print the identifier (if necessary)
# print the upper/lower identifier, if set
    if ( defuplo != "" ) {
	print defuplo "_BOUND_DATA"
	defuplo = ""
    }
# take deftyp value if no typ was found
    if ( typ == "none" ) typ = deftyp
# default initial value set to NOE
    if ( deftyp == "none" ) {
	if ( typ == deftyp ) typ = "NOE"
	deftyp = typ
    }
# check against type last determined from a comment or identifier (if any)
    if ( typ != deftyp && deffrom != "none" ) {
	print "* Conflict: last", deffrom, "says <" deftyp ">, but data interpreted as <" typ ">" >> logfil
	print "* line", NR ":", org >> logfil
	if ( fatal ) exit( 1 )
# this is to prevent too many errors
	if ( !nosavtyp ) {
	    deftyp = typ
	    deffrom = "data"
	}
    }
# print the derived identifier
    if ( typ != dattyp ) {
	print "START_OF_" typ "_DATA"
	dattyp = typ
    }
}

function prtsto() {
    for ( i=1; i<=cmmcnt; i++ )
	print cmmsto[i]
    cmmcnt = 0
}

function prtlin() {
    prtsto()
    if ( dismanreslin != "" ) {
	print dismanreslin
# the line with the first residue number and name will be repeated every time
# but in this way the switch between NOE and HB at least works
#	dismanreslin = ""
    }
    print org
}

function prtend() {
    prtsto()
    if ( ide ) print "END_OF_DATA"
}

function prtmodel() {
# write a model identifier only if the previous one occurred more 
# than 10 lines ago, and exit at second model
    if ( modelline < ( NR - 10 ) ) {
	if ( ++modelnum > 1 ) exit
##	if ( ide ) print "START_OF_MODEL"
	print "START_OF_MODEL"
    }
    modelline = NR
}


# START PROCESSING - skip blank lines - change to uppercase
# reset the type flag and atom names used in readatXYZ and guessdat

NF == 0 { 
    cmmsto[++cmmcnt] = $0
##    prtsto()
    next
}

{
    org = $0
    $0 = toupper( $0 )
    AA = ""
    BB = ""
    typ = "none"
}

# Start of coordinate section of "mr" file

mr  {
    if ( $0 ~ /^ *START OF MODEL/ || $0 ~ /^ *ENTRY H ATOM NAME/ || \
         $0 ~ /^ *RAW FILE HAD/ )  {
	if ( !modelnum ) prtend()
	prtmodel()
    }
}

# Now writing model

modelnum {
    print org
    next
}

# Preset identifiers (may be preceded by comment delimiter ! or #)
# keep ide value (comment/data interpretation may still override)
# unset cmm in unspecified data blocks

/^[#!* ]*END_OF_DATA *$/ {
    cmm = 0
    sub( /^[#! ]*/, "", org )
    prtlin()
    next
}

/^[#!* ]*START_OF_.*DATA *$/ {
    sub( /^[#! ]*START_OF_/, "" )
    sub( /_*DATA *$/, "" )
    if ( $0 == "" ) $0 = "OTH"
    if ( $0 == "OTH" )
	cmm = 0
    else
	cmm = cmmorig
    if ( ide ) {
	deftyp = $0
	deffrom = "identifier"
	dattyp = $0
    }
    sub( /^[#! ]*/, "", org )
    prtlin()
    next
}

/^[#!* ]*START_OF_MODEL *$/ {
# do exactly the same as above (Start of coordinate section of "mr" file)
    if ( !modelnum ) prtend()
    prtmodel()
    if ( !ide ) {
	sub( /^[#! ]*/, "", org )
	prtlin()
    }
    next
}

# Comment lines

/^ *[#!*]/ {
    guesscmm()
    cmmsto[++cmmcnt] = org
    next
}

## JFD 10/04/2001 adds to deal with entry 1eii special case code
## This was causing the DIANA script to wrongly convert some lines into
## restraints and that caused the count to be wrong in the .noe file
## tracked by my Python scripties.
/^CONECT/ {
    garbage()
    next
}

# Now try to determine the type of data by text interpretation - EVERY line

# If the data are "OTH", print them immediately by calling prtide and prtlin, 
# just as for UNKNOWN file type (see below)
{
    guesscmm()
    if ( deftyp == "OTH" ) {
	prtide()
	prtlin()
	next
    }
}

# Process other data types depending on the file type
# Functions to call: readatXYZ, guessdat (calling prtide), prtlin
#                or: prtide, prtlin
 
# DISGEO type (special case: identifiers not needed)

filtyp == 1 && ( \
/^NOEUPP/ || /^HBUPP/  || /^SSUPP/ || \
/^NOELOW/ || /^HBLOW/  || /^SSLOW/ || \
/^ANGLE/  || /^JCOUPL/ || /^SHIFTS/ ) { 
    prtlin()
    next
}

# XPLOR type

filtyp == 2 {
    if ( $0 ~ /^ *ASSI/ || $0 ~ /RESI.*AND.*NAME/ ) {
	lastdata = NR
	readatXPLOR()
	guessdat()
	prtlin()
	next
    }
    if ( $0 ~ /^ *REMARK/ || $0 ~ /^ *EVALUATE/ || $0 ~ /^ *SET / || $0 ~ /^ *END *$/ ) {
	garbage()
	next
    }
# try to filter out "broken data lines"
    if ( NR - lastdata <= 2 ) {
	print org
	next
    }
}

# BIOSYM type

filtyp == 3 {
# strictly speaking non-NOE type restraints can only occur as "general distance restraints"
#    if ( deftyp == "DIS" )
	readatBIOSYM()
    guessdat()
    prtlin()
    next
}

# DISMAN/DIANA type

filtyp == 4 {    
    nfield = countdatafields()
    if ( deftyp == "TOR" ) {
	if ( (nfield==5||nfield==6) && $1 ~ /^[0-9]+$/ ) {
	    dismanreslin = ""
	    prtide()
	    prtlin()
	    next
	}
    }
    else {
	if ( ( nfield == 2 && $1 ~ /^[0-9]+$/ ) || ( (nfield==5||nfield==6) && $2 ~ /^[0-9]+$/ && $5 ~ /^[0-9.]+$/ ) || ( (nfield==7||nfield==8) && $1 ~ /^[0-9]+$/ && $4 ~ /^[0-9]+$/ && $7 ~ /^[0-9.]+$/ ) ) {
	    if ( nfield == 2 ) {
		dismanreslin = org
		next
	    }
	    if ( nfield > 6 ) {
		dismanreslin = ""
	    }
	    readatDISMAN()
	    guessdat()
	    prtlin()
	    next
	}
    }
# non-fitting lines will be handled by garbage
# it's not easy to filter them out, because of dismanreslin
}

# UNKNOWN type 

filtyp == filtyp_unk {
    guessdatUNK()
    prtlin()
    next
}

# All other lines: treat as comment

    {
    garbage()
}

END {
    prtsto()
}
